import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm

# 原始数据
y = np.array([
    11.2, 13.4, 40.7, 5.3, 24.8, 12.7, 20.9, 35.7, 8.7, 9.6,
    14.5, 26.9, 15.7, 36.2, 18.1, 28.9, 14.9, 25.8, 21.7, 25.7,
])

X = np.array([
    [587000, 16.5, 6.2],
    [643000, 20.5, 6.4],
    [635000, 26.3, 9.3],
    [692000, 16.5, 5.3],
    [1248000, 19.2, 7.3],
    [643000, 16.5, 5.9],
    [1964000, 20.2, 6.4],
    [1531000, 21.3, 7.6],
    [713000, 17.2, 4.9],
    [749000, 14.3, 6.4],
    [7895000, 18.1, 6.0],
    [762000, 23.1, 7.4],
    [2793000, 19.1, 5.8],
    [741000, 24.7, 8.6],
    [625000, 18.6, 6.5],
    [854000, 24.9, 8.3],
    [716000, 17.9, 6.7],
    [921000, 22.4, 8.6],
    [595000, 20.2, 8.4],
    [3353000, 16.9, 6.7],
])

# 添加截距项
X_with_intercept = sm.add_constant(X)

# 使用 statsmodels 进行回归分析
model = sm.OLS(y, X_with_intercept).fit()

print("多元线性回归结果:")
print(f"截距 (Intercept): {model.params[0]:.6f}")
print("\n估计系数 (Coefficients):")
var_names = ["Inhabitants", "Percent with incomes below $5000", "Percent unemployed"]
for i, var_name in enumerate(var_names):
    print(f"  {var_name}: {model.params[i+1]:.6f}")

print("\n标准误差 (Standard Errors):")
for i, var_name in enumerate(var_names):
    print(f"  {var_name}: {model.bse[i+1]:.6f}")

print("\nt 值 (t-values):")
for i, var_name in enumerate(var_names):
    print(f"  {var_name}: {model.tvalues[i+1]:.6f}")

# 计算 VIF
print("\n方差膨胀系数 (VIFs):")
vif_values = []
for i in range(X.shape[1]):
    vif = variance_inflation_factor(X_with_intercept, i+1)  # 跳过常数项
    vif_values.append(vif)
    print(f"  {var_names[i]}: {vif:.6f}")

print(f"\nR²: {model.rsquared:.6f}")
print(f"调整 R²: {model.rsquared_adj:.6f}")

print("\nVIF 解释:")
print("  VIF < 5: 无显著多重共线性")
print("  5 ≤ VIF < 10: 中等多重共线性") 
print("  VIF ≥ 10: 严重多重共线性")
for i, vif in enumerate(vif_values):
    level = "无显著多重共线性"
    if vif >= 10:
        level = "严重多重共线性"
    elif vif >= 5:
        level = "中等多重共线性"
    print(f"  {var_names[i]}: {level} (VIF={vif:.2f})")

print("\n完整回归摘要:")
print(model.summary())